Intro to Bootstrap

Edward Vytlacil

Application: Mincer Wage Equation

  • Application: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Average Effect of Experience

  • Experience level that maximizes log wages.

  • Bootstrap Hypothesis Test

  • Other Bootstrap Variants

Application: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

  • For ease of exposition, dropping education.

Application: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

  • Average effect of experience: \[\beta_1 + 2 \cdot \beta_2 \cdot \mathbb{E}(\mbox{exp}),\]

  • Experience level that maximizes log wage: \[\left | \frac{\beta_1}{2 \beta_2}\right |.\]

Application: Mincer Wage Equation

  • Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

  • Let \(\overline{\mbox{exp}}_N\) denote sample mean of experience.

  • Estimate average effect of experience by \[\widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N,\]

    • Estimator is a non-linear function of \((\widehat \beta_{1,N}, \widehat \beta_{2,N}, \overline{\mbox{exp}}_N)\).

Application: Mincer Wage Equation

  • Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

  • Estimate experience level that maximizes log wage by \[\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right|.\]

    • Estimator is a non-linear functions of \((\widehat \beta_{1,N}, \widehat \beta_{2,N})\).

Application Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
  • Consider estimation using CPS 1985 data.

    • available in package AER.

    • after loading AER, see ?CPS1985 for description of data.

  • Consider only using observations with \(12\) years of schooling.

Application Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
[1] "Mean Experience: 18.14"

Call:
lm(formula = I(log(wage)) ~ experience + I(experience^2), data = df)

Coefficients:
    (Intercept)       experience  I(experience^2)  
      1.5829802        0.0352680       -0.0005831  

Application Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp
[1] "Estimated Avg Effect Experience: 0.0141"
[1] "Estimated Experience that Maximizes Wage: 30.2"

Application Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp
  • \(\overline{\mbox{exp}}_N=\) 18.1

  • \(\widehat \beta_{1,N} =\) 0.0353

  • \(\widehat \beta_{2,N} =\) -0.000583

  • \(\widehat \beta_{1,N} + 2 \widehat \beta_{2,N} \overline{\mbox{exp}}_N=\) 0.0141

  • \(\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right| =\) 30.2

Application: Mincer Wage Equation

  • \(\widehat \beta_{1,N} + 2 \widehat \beta_{2,N} \overline{\mbox{exp}}_N=\) 0.0141

  • \(\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right| =\) 30.2

  • How to compute standard errors? construct confidence intervals? perform hypothesis tests?

Delta Method vs Bootstrap

  • Application: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Average Effect of Experience

  • Experience level that maximizes log wages.

  • Bootstrap Hypothesis Test

  • Other Bootstrap Variants

Traditional approach: Delta Method.

  • Asymptotic approximation based on linear approximation from Mean Value Theorem combined with LLN, CLT.

  • Limitations:

    • only justified asymptotically,
    • sometimes provides poor approximation,
    • derivations using delta method can be tedious, chance of algebraic errors.

Bootstrap principle

  • Treat sample of size \(N\) as if it were the population.

  • Perform MC simulation, creating many bootstrap samples of size \(N\) by resampling with replacement from original sample.

  • Compute estimator on each sample, use distribution of estimates across samples to approximate distribution of estimator.

  • Difference between drawing from population vs from sample should vanish as \(N \rightarrow \infty\).

Bootstrap: Basic Principle

Bootstrap Method

  • Limitations:

    • only justified asymptotically,
    • sometimes provide poor approximation,
    • can be computationally intensive
  • Advantages:

    • Can bypass derivations, complicated formulas. Less tedious to implement. Removes chance of algebraic error in derivations.
    • Sometimes provides better approximation than delta method.

Bootstrap Method

  • There are many versions of the bootstrap, but we will focus on nonparametric bootstrap and corresponding . . .

    • estimate of bias,
    • standard errors,
    • normal-approximation CI,
    • percentile CI,
    • percentile-t CI,
    • hypothesis tests.

Bootstrap Algorithm

  • Application: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Average Effect of Experience

  • Experience level that maximizes log wages.

  • Bootstrap Hypothesis Test

  • Other Bootstrap Variants

Sample and Estimator

  • Suppose \(\mathbf{W} = (W_1,W_2,...,W_N)\) is observed sample
    • e.g., \(W_i = (Y_i,X_i)\) in regression context.
  • Estimator \(\widehat \theta = s(\mathbf{W})\)
    • expressing estimator as function of data, rule for how to use data to compute estimate.
    • in our application, \(\widehat \theta\) is estimator of average effect of experience or of experience level that maximizes wages.

Bootstrap algorithm

Bootstrap Samples

Bootstrap Replications

  • Create \(B\) bootstrap samples, each constructed by drawing from original sample with replacement \(N\) times.

  • Compute bootstrap replication on each bootstrap sample, \(\widehat \theta^*(b) = s(\mathbf{W}^{*b})\)

\[\begin{align} \mathbf{W}^{*1} & = (W^{*1}_1,W^{*1}_2,…,W^{*1}_N) ~\Rightarrow \widehat \theta^*(1)=s(\mathbf{W}^{*1}) \\ \mathbf{W}^{*2} & = (W^{*2}_1,W^{*2}_2,…,W^{*2}_N) ~ \Rightarrow \widehat \theta^*(2)=s(\mathbf{W}^{*2}) \\ & \qquad \qquad \quad ~~~\vdots \qquad \qquad \qquad ~~~ \vdots \\ \mathbf{W}^{*B} & = (W^{*B}_1,W^{*B}_2,…,W^{*B}_N) \Rightarrow \widehat \theta^*(B)=s(\mathbf{W}^{*B}). \\ \end{align}\]

  • Use empirical distribution of \((\widehat \theta^*(1), \widehat \theta^*(2),...,\widehat\theta^*(B))\) to approximate distribution of \(\widehat \theta\).

Example: Algorithm for Bootstrap s.e.

  • Create \(B\) bootstrap samples, each constructed by drawing from original sample with replacement \(N\) times.

  • Compute bootstrap replication \(\widehat \theta^*(b)\) on each bootstrap sample.

  • Define bootstrap s.e. for \(\widehat \theta\) by \[\widehat{\mbox{se}}_B=\mbox{s.d.}(\widehat \theta^*(1),\widehat \theta^*(2),…,\widehat\theta^*(B))\] (illustrated in figure to left).

Normal Approx. Bootstrap CI

Given bootstrap standard error, \(\widehat{\mbox{se}}_B\):

  • Construct 95% normal-approximation bootstrap CI on \(\theta\) by \[C^{nb} = [\widehat{\theta}-1.96 \cdot \widehat{\mbox{se}}_B, ~\widehat{\theta}+1.96 \cdot \widehat{\mbox{se}}_B]\]

  • More generally, construct \(1-\alpha\) level normal-approximation bootstrap CI on \(\theta\) by \[[\widehat{\theta}- z_{1-\alpha/2} \cdot \widehat{\mbox{se}}_B, ~\widehat{\theta}+z_{1-\alpha/2} \cdot \widehat{\mbox{se}}_B]\] where \(z_{1-\alpha/2}\) is \(1-\alpha/2\) quantile of standard normal distribution.

Example: Bootstrap Estimate Bias

  • Given \(B\) bootstrap replications, define the boostrap estimate of bias by \[\widehat{\mbox{bias}}^*= \frac{1}{B} \sum_{b=1}^B \widehat \theta^*(b) - \widehat \theta \] where \(\widehat \theta\) is estimate from original sample.

  • why centered at \(\widehat \theta\)?

Example: Algorithm for Percentile CI

  • Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((\widehat \theta^*(1),\widehat \theta^*(2),…,\widehat\theta^*(B))\)

    • compute \(q^*_{\alpha}\) as the value such that fraction \(\alpha\) of the bootstrap replications are less than \(q^*_{\alpha}\).
  • Bootstrap Percentile CI: \[C^{PC}= [q_{\alpha/2},~ q_{1-\alpha/2}].\]

Example: Algorithm for Percentile-t CI

  • Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((T^*(1), T^*(2), ...,T^*(B))\), where \[T^*(b) = \frac{\widehat \theta^*(b) - \widehat \theta}{s(\widehat \theta^*(b))}\]

where \(s(\widehat \theta^*(b))\) is the standard error of \(\widehat \theta^*(b)\) computed using only the \(b\)th bootstrap replication.

  • Note centering at \(\widehat \theta\), and note s.e. using only \(b\)th bootstrap replication.

Example: Algorithm for Percentile-t CI (cont’d)

  • Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((T^*(1), T^*(2), ...,T^*(B))\).

  • Bootstrap Percentile-t CI: \[C^{Pt}= [\widehat \theta - s(\widehat \theta) q_{1-\alpha/2},~ \widehat \theta - s(\widehat \theta) q_{\alpha/2}].\]

Pros and Cons Bootstrap Methods

  • Bootstrap standard errors:
    • Pros:
      • intuitive, easy to implement,
      • researchers familiar with s.e., is the norm to report s.e..
    • Cons:
      • Valid only under relatively strong conditions,
      • May be erratic, poorly behaved.

Pros and Cons Bootstrap Methods

  • Normal-approximation bootstrap CI:

    • Pros:
      • Intuitive, easy to implement,
      • Researchers familiar with normal CIs,
      • Many estimators are asymptotically normal.

Pros and Cons Bootstrap Methods

  • Normal-approximation bootstrap CI:

    • Cons:
      • Invalid if bootstrap s.e. are invalid,
      • Poorly behaved if bootstrap s.e. are poorly behaved,
      • Normality may be bad approximation
        • especially poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

  • bootstrap percentile CI:

    • Pros:

      • intuitive, easy to implement,

      • valid under weaker conditions than bootstrap s.e. or bootstrap normal CI.

    • Cons:

      • can provide poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

  • bootstrap percentile-t CI:

    • Pros:

      • valid under weaker conditions than bootstrap s.e. or bootstrap normal CI.

      • better theoretical properties (more accurate) than percentile CI.

    • Cons:

      • need to compute standard error on each bootstrap replication.

Average Effect of Experience

  • Application: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Average Effect of Experience

  • Experience level that maximizes log wages.

  • Bootstrap Hypothesis Test

  • Other Bootstrap Variants

Avg Effect of Experience

  • We now implement bootstrap s.e., CIs and hypothesis tests for average effect of experience.

    • \(\theta = \beta_{1} + 2 \cdot \beta_{2} \cdot \mathbb{E}[\mbox{exp}],\)

    • \(\widehat \theta = \widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N.\)

  • We first proceed in base R, then consider using package boot.

  • First step: create bootstrap samples.

Avg Effect of Experience

  • Create bootstrap sample:

    • Randomly draw \(N\) integers, with replacement, from \((1,...,N)\), with each integer equally likely to be draw,

    • Create new sample formed by observations with those indices.

    • Some observations from original sample will be drawn multiple times into bootstrap sample, others will not appear in bootstrap sample.

N <- nrow(df) 
indices.b1<-sample(1:N,replace=TRUE)  # draw indices with replacement from (1,...,N), N sample size
tab.b1 <- df[indices.b1,] # construct bootstrap sample

Bootstrap Sample 1

Original Sample

Bootstrap Sample 1

Bootstrap Sample 2

Original Sample

Bootstrap Sample 2

Bootstrap Sample 3

Original Sample

Bootstrap Sample 3

Bootstrap Sample 4

Original Sample

Bootstrap Sample 4

Bootstrap Replications Avg. Effect Experience

  • given bootstrap sample, \(\mathbf{W}^{*b}=(W^{*b}_1,W^{*b}_2,…,W^{*b}_N)\), compute \(\widehat \theta^*(b)=s(\mathbf{W}^{*b}).\)

    • compute mean experience on bootstrap sample, \(\overline{\mbox{exp}}^{~*b}\).
    • run regression on bootstrap sample, obtain bootstrap sample OLS coefficients \((\widehat \beta_{0}^{*b},\widehat \beta_{1}^{*b}, \widehat \beta_{2}^{*b})\)
    • compute \[\widehat \theta^*(b) = \widehat \beta_{1}^{*b} + 2 \cdot \widehat \beta_{2}^{*b} \cdot \overline{\mbox{exp}}^{~*b}.\]

Bootstrap Replications Avg. Effect Experience

#  tab.b1 <- df[indices.b1,], bootstrap sample
mean.exp.b1 <- mean(tab.b1$experience) #  sample mean experience on bootstrap sample
reg.b1 <- lm(I(log(wage))~ experience+I(experience^2),data=tab.b1) # OLS on bootstrap sample
# form bootstrap replication of estimated mean effect experience
a <- c(0,1,2*mean.exp.b1)
theta.b1 <- as.numeric(t(a)%*%coef(reg.b1)) 

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 1, \(\widehat{\theta}^*(1) =\) 0.032 \(-2\cdot\) -0.0005 \(\cdot\) 18.05

\(=\) 0.0139

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
18.05

\(\widehat{y}=\) 1.64 \(+\) 0.032 \(\mbox{exp}\) \(-\) 0.0005 \(\mbox{exp}^2\)

\(\widehat{\theta} =\) 0.035 \(+2\cdot\) 0.0353 \(\cdot\) 18.14 \(=\) 0.0141

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 2, \(\widehat{\theta}^*(2) =\) 0.024 \(-2\cdot\) -0.00044 \(\cdot\) 17.17

\(=\) 0.0085

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
17.17

\(\widehat{y}=\) 1.72 \(+\) 0.024 \(\mbox{exp}\) \(-\) 0.00044 \(\mbox{exp}^2\)

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 3, \(\widehat{\theta}^*(3) =\) 0.023 \(-2\cdot\) -0.00034 \(\cdot\) 17.89

\(=\) 0.0107

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
17.89

\(\widehat{y}=\) 1.64 \(+\) 0.023 \(\mbox{exp}\) \(-\) 0.00034 \(\mbox{exp}^2\)

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 4, \(\widehat{\theta}^*(4) =\) 0.046 \(-2\cdot\) -0.00078 \(\cdot\) 19.19

\(=\) 0.0159

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
19.19

\(\widehat{y}=\) 1.5 \(+\) 0.046 \(\mbox{exp}\) \(-\) 0.00078 \(\mbox{exp}^2\)

Construct \(B\) Bootstrap Replications

  • We have created 4 bootstrap samples, and corresponding 4 bootstrap replications.

  • We want to create \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications, for \(B\) large, for example, \(B=10,000\).

  • How to repeat steps many times?

    • we first consider approach in base R using replicate,
    • we then consider using functions from package boot.

Use replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))
  • Create user defined function and use replicate to generate \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))
  • User defined function takes two arguments:
    1. original sample as data frame,
    2. indices;
  • Function returns bootstrap replication.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))
  • Inputting user defined function, original data and random sample of indices into replicate function, which is returning \(B\) bootstrap replications.

  • Use output to compute bootstrap s.e. or CI.

Workflow using base R

flowchart LR
  A{     User-Defined Fn.    <br>     for Estimator    <br>   } --> B{       `replicate`  }
  E{     Original     <br>    Sample    <br>   } --> B 
 B --> C{    bootstrap s.e    <br>}
 B --> D{  bootstrap CI    <br>}

Distribution of Bootstrap Estimates

Here: \[\begin{align} \theta = ~ & \beta_{1} + 2 \cdot \beta_{2} \cdot \mathbb{E}(\mbox{exp})\\ \widehat \theta = ~ & \widehat \beta_{1} + 2 \cdot \widehat \beta_{2} \cdot \overline{\mbox{exp}}\\ \widehat \theta^*(b) = ~ & \widehat \beta_{1}^{*b} + 2 \cdot \widehat \beta_{2}^{*b} \cdot \overline{\mbox{exp}}^{~*b} \end{align}\]

Key Idea:

  • Use simulated distribution of bootstrap sample statistics \(\widehat \theta^*(b)\) to approximate distribution of \(\widehat \theta\).

Distribution of Bootstrap Estimates

Graphing histogram of \((\widehat \theta^*(1), \widehat \theta^*(2), ... , \widehat \theta^*(B))\).

fig.boot <- ggplot(data.frame(boot.est), 
            aes(x=boot.est))  +
    xlab("Bootstrap Replications: Avgerage Effect Experience")+
    geom_histogram(color="darkblue", 
            fill="lightblue",binwidth=0.0005)
fig.boot 

Key Idea:

  • Use simulated distribution of bootstrap sample statistics \(\widehat \theta^*(b)\) to approximate distribution of \(\widehat \theta\).

Bootstrap Estimate Bias

Bootstrap estimate bias:

\(\bar{\theta}^* = \frac{1}{B} \sum_{b=1}^B \widehat \theta^{}(b)\)

\(\widehat{\mbox{bias}}^* = \bar{\theta}^* - \widehat \theta\)

  • Use difference between mean of bootstrap replications and estimate on original sample as bootstrap estimator of bias.

Bootstrap Estimate Bias

Bootstrap estimate bias:

mean(boot.est)-theta.hat
[1] 0.000051192
  • Use difference between mean of bootstrap replications and estimate on original sample as bootstrap estimator of bias.

    • \(\widehat{\mbox{bias}}^*=\) 0.0000512

Bootstrap Standard Errors

\[\mbox{s.e.}^{\mbox{boot}}(\widehat{\theta}_N)=\sqrt{ \frac{1}{B-1} \sum_{b=1}^B (\widehat \theta^{*}(b) - \bar{\theta}^*)^2 }\]

sd(boot.est)
[1] 0.002900513
  • Use s.d. of bootstrap replications \(\widehat \theta^*(b)= \widehat \beta_{1}^{*}(b) + 2 \cdot \widehat \beta_{2}^*(b) \cdot \overline{\mbox{exp}}^*(b)\) to form bootstrap standard errors of \(\widehat \theta = \widehat \beta_{1} + 2 \cdot \widehat \beta_{2} \cdot \overline{\mbox{exp}}\).

    • \(\mbox{s.e.}^{\mbox{boot}}(\widehat \theta)=\) 0.0029

Normal Approx. Bootstrap CI

95% normal-approximation bootstrap CI: \[C^{nb} = [\widehat{\theta}-1.96 \cdot \widehat{\mbox{se}}_B, ~\widehat{\theta}+1.96 \cdot \widehat{\mbox{se}}_B]\] Based on approximate normality.

  • Plug bootstrap s.e. into conventional 95% CI for asymptotically normal estimators.

Normal Approx. Bootstrap CI

Endpoints of Normal Bootstrap CI:

theta.hat -qnorm(.975)*sqrt(var(boot.est))
[1] 0.008424728
theta.hat +qnorm(.975)*sqrt(var(boot.est))
[1] 0.01979453
  • Plug bootstrap s.e. into conventional 95% CI for asymptotically normal estimators:
    • \(C^{nb} =\) [ 0.00842 \(,\) 0.0198 ].

Bootstrap Percentile CI

95% Bootstrap Percentile Interval: \[ C^{pc} = [q^*_{0.025}, ~q^*_{0.975}].\] where \(q^*_{0.025}\) and \(q^*_{0.975}\) are the \(0.025\) and \(0.975\) empirical quantiles of bootstrap replications.

Bootstrap Percentile CI

Bootstrap quantiles:

q <- quantile(boot.est, c(0.025, 0.975)) 
q 
      2.5%      97.5% 
0.00852697 0.01994357 
  • Bootstrap Percentile 95% CI:
    \(C^{pc} =\) [ 0.00853, 0.0199 ]

Boot. Percentile vs Norm Approx CI

  • Bootstrap Normal Approx. 95% CI:

    • [ 0.00842 , 0.0198 ].
  • Bootstrap Percentile 95% CI:

    • [ 0.00853 , 0.0199 ].
  • Why so similar?

Alternatively Can Use Package boot

Workflow using package boot.

flowchart LR
  A{      User-Defined Fn.    <br>    for Estimator   } --> B{        `boot`   }
  F{     Original Sample    } --> B 
 B --> C{    bootstrap s.e    }
  B --> D{    `boot.ci`   }
  D --> E{    bootstrap CI    }

Alternatively using package boot

Function boot from package boot:

parameter description
data original sample (as vector, matrix, or data frame)
statistic function that constructs bootstrap replication on bootstrap sample, must have two arguments: (1) data, and (2) indices.
R number of bootstrap samples/replications, we have been denoting \(B\)
stype for this course, will always set to "i" for indices.

Alternatively using package boot

Function boot.ci from package boot:

parameter

description

boot.out

output from boot

conf

confidence level, by default is 0.95.

type

  • “norm” for normal approx. CI,

  • “perc” for percentile CI,

  • “bca” for bootstrap accelerated bias-corrected percentile interval.

  • stud for bootstrap percentile-t CI

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • user defined function takes two inputs:
    • data frame (original sample)
    • indices;
  • function returns bootstrap replication.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • boot taking as inputs:

    • original sample (data frame),
    • user defined function,
    • number of replications.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • to understand output from boot, enter:

    • ?boot
    • str(boot.out1)

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot1, R = B, stype = "i")


Bootstrap Statistics :
      original        bias    std. error
t1* 0.01410963 0.00004975603 0.002904684

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • boot.ci taking as inputs:

    • saved output from boot,
    • choice of confidence level,
    • type of bootstrap CI.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • to understand output from boot.ci, enter:

    • ?boot.ci
    • str(boot.ci1)

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out1, conf = 0.95, type = c("norm", "perc", 
    "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   ( 0.0084,  0.0198 )   ( 0.0085,  0.0199 )   ( 0.0083,  0.0197 )  
Calculations and Intervals on Original Scale

Output from function boot

boot.out1<- boot(data=df, statistic=f.boot1, 
                 R=B, stype="i")

#boot.out1$t0 is estimate on original sample
boot.out1$t0
[1] 0.01410963
#boot.out1$t is vector of 10,000 stored 
#            bootstrap replications

ggplot(data.frame(boot.out1$t), 
            aes(x=boot.out1$t))  +
    xlab("Bootstrap Replications: 
            Avgerage Effect Experience")+
    geom_histogram(color="darkblue", 
            fill="lightblue",binwidth=0.0005)

plot(boot.out1)

summary(boot.out1$t)
       V1          
 Min.   :0.003135  
 1st Qu.:0.012088  
 Median :0.014050  
 Mean   :0.014060  
 3rd Qu.:0.015995  
 Max.   :0.024635  
summary(boot.out1)
      R original     bootBias    bootSE bootMed
1 10000  0.01411 -0.000049418 0.0028694 0.01405

Analytical s.e.

Let \(\sigma^2_X = \mbox{Var}(X_i)\), and

\[ a = ~ \left[ \begin{array}{c} 0\\ 1 \\ 2 \mathbb{E}(\mbox{exp}) \end{array} \right]^{\prime}, \qquad\qquad \widehat a= \left[ \begin{array}{c} 0\\ 1 \\ 2 \overline{\mbox{exp}} \end{array} \right]^{\prime}. \]

Analytical s.e. ignoring sampling variability in \(\overline{\mbox{exp}}\), \[\mbox{s.e.}(\widehat \theta) = \frac{1}{\sqrt{N}} \left( \widehat a^{\prime} \widehat \Sigma \widehat a \right)^{1/2}\]

Analytical s.e.

Let \(\sigma^2_X = \mbox{Var}(X_i)\), and \[ a = ~ \left[ \begin{array}{c} 0\\ 1 \\ 2 \mathbb{E}(\mbox{exp}) \end{array} \right]^{\prime}, \qquad\qquad \widehat a= \left[ \begin{array}{c} 0\\ 1 \\ 2 \overline{\mbox{exp}} \end{array} \right]^{\prime}. \]

Analytical s.e. using delta method, \[\mbox{s.e.}(\widehat \theta) = \frac{1}{\sqrt{N}} \left( \widehat a^{\prime} \widehat \Sigma \widehat a + 4 \widehat \beta_2^2 \widehat \sigma^2_X \right)^{1/2}.\]

Analytical s.e.

# not accounting for variability in sample mean
library(sandwich)
a <- c(0,1,2*mean.exp)
se.1 <-  as.numeric(sqrt(t(a)%*%vcovHC(reg.1)%*%a))
signif(se.1,3)
[1] 0.00277
# corresponding 95% CI
c(signif(theta.hat -qnorm(.975)*se.1,3), signif(theta.hat+qnorm(.975)*se.1,3))
[1] 0.00869 0.01950
# using delta method
se.2 <- as.numeric(sqrt(t(a)%*%vcovHC(reg.1)%*%a +
                          4*(coef(reg.1)[[3]])^2*var(df$experience)/nrow(df)))
signif(se.2,3)
[1] 0.00293
# corresponding 95% CI
c(signif(theta.hat -qnorm(.975)*se.2,3), signif(theta.hat +qnorm(.975)*se.2,3))
[1] 0.00836 0.01990

Bootstrap Percentile-t CI

  • For percentile-t CI, need to estimate variance of estimator on each bootstrap sample.

  • Can do so using result from Delta method.

Bootstrap Percentile-t CI

f.boot1.t <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  est.b <- as.numeric(t(a)%*%coef(reg.b)) 
  # form estimated variance
  var.b <- as.numeric(t(a)%*%vcovHC(reg.b)%*%a+4*(coef(reg.b)[[3]])^2*var(d[i,]$experience)/nrow(d[i,]) )
  c(est.b,var.b)
  }
boot.out1.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")
  • Notice that function f.boot1.t is returning variance of estimator computed on bootstrap sample as second component, necessary if inputting into boot for later input into boot.ci to construct percentile-t CI.

Bootstrap Percentile-t CI

f.boot1.t <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  est.b <- as.numeric(t(a)%*%coef(reg.b)) 
  # form estimated variance
  var.b <- as.numeric(t(a)%*%vcovHC(reg.b)%*%a+4*(coef(reg.b)[[3]])^2*var(d[i,]$experience)/nrow(d[i,]) )
  c(est.b,var.b)
  }
boot.out1.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")
summary(boot.out1.t$t)
       V1                 V2             
 Min.   :0.002754   Min.   :0.000004574  
 1st Qu.:0.012117   1st Qu.:0.000007711  
 Median :0.014067   Median :0.000008607  
 Mean   :0.014086   Mean   :0.000008716  
 3rd Qu.:0.016044   3rd Qu.:0.000009618  
 Max.   :0.025135   Max.   :0.000015271  

V2 is estimated variance on bootstrap sample.

Bootstrap Percentile-t CI

boot.out1.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")
boot.ci1.t <- boot.ci(boot.out1.t, conf=0.95, type=c("norm","perc","bca","stud"))
boot.ci1.t
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out1.t, conf = 0.95, type = c("norm", 
    "perc", "bca", "stud"))

Intervals : 
Level      Normal            Studentized     
95%   ( 0.0083,  0.0199 )   ( 0.0084,  0.0199 )  

Level     Percentile            BCa          
95%   ( 0.0085,  0.0199 )   ( 0.0084,  0.0199 )  
Calculations and Intervals on Original Scale

Summary Avg Effect Exp

Estimated Bias -0.0000494
Bootstrap s.e. 0.00287
Analytical s.e. (Delta method) 0.00293
Bootstrap Normal 95% CI [0.00837, 0.0198]
Normal 95% CI (Delta method) [0.00836, 0.0199]
Bootstrap Percentile 95% CI [0.00849, 0.0199]
BCa 95% CI [0.0083, 0.0197]
Bootstrap Percentile-t 95% CI [0.00839, 0.0199]

How to summarize bootstrap results?

What results would you report? how to interpret?

Rerun Analysis?

  • Can rerun bootstrap simulation to see if results are stable across simulations.
  • If not?
    • insufficient number of replications? increase \(B\)?
    • problem with bootstrap method in particular example?
      • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.

Rerun Analysis?

  • Can rerun bootstrap simulation to see if results are stable across simulations.
  • If not?
    • problem with bootstrap method in particular example?
      • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.
      • percentile CIs are valid under weaker assumptions and tend to be more robust/stable.

Repeated Simulations

boot.out2.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")
plot(boot.out2.t)

summary(boot.out2.t)[1,]
      R original    bootBias    bootSE  bootMed
1 10000  0.01411 0.000021084 0.0029555 0.014118
summary(boot.out2.t$t)
       V1                 V2             
 Min.   :0.003235   Min.   :0.000004581  
 1st Qu.:0.012130   1st Qu.:0.000007707  
 Median :0.014118   Median :0.000008596  
 Mean   :0.014131   Mean   :0.000008727  
 3rd Qu.:0.016142   3rd Qu.:0.000009569  
 Max.   :0.026139   Max.   :0.000016309  
boot.out3.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")
plot(boot.out3.t)

summary(boot.out3.t)[1,]
      R original     bootBias    bootSE  bootMed
1 10000  0.01411 -0.000034136 0.0029183 0.014063
summary(boot.out3.t$t)
       V1                 V2             
 Min.   :0.002633   Min.   :0.000004569  
 1st Qu.:0.012082   1st Qu.:0.000007673  
 Median :0.014063   Median :0.000008591  
 Mean   :0.014075   Mean   :0.000008711  
 3rd Qu.:0.015980   3rd Qu.:0.000009593  
 Max.   :0.026173   Max.   :0.000015257  

Repeated Simulations

boot.ci(boot.out2.t, conf=0.95, type=c("norm","perc","bca","stud")) #2nd Simulation 
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out2.t, conf = 0.95, type = c("norm", 
    "perc", "bca", "stud"))

Intervals : 
Level      Normal            Studentized     
95%   ( 0.0083,  0.0199 )   ( 0.0085,  0.0199 )  

Level     Percentile            BCa          
95%   ( 0.0084,  0.0199 )   ( 0.0083,  0.0199 )  
Calculations and Intervals on Original Scale
boot.ci(boot.out3.t, conf=0.95, type=c("norm","perc","bca","stud")) # 3rd Simulation
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out3.t, conf = 0.95, type = c("norm", 
    "perc", "bca", "stud"))

Intervals : 
Level      Normal            Studentized     
95%   ( 0.0084,  0.0199 )   ( 0.0085,  0.0199 )  

Level     Percentile            BCa          
95%   ( 0.0085,  0.0198 )   ( 0.0085,  0.0199 )  
Calculations and Intervals on Original Scale

Summary Repeated Simulations

Simulation 1 Simulation 2 Simulation 3
Estimated bias -0.0000494 0.0000211 -0.0000341
Bootstrap s.e. 0.00294 0.00296 0.00292
Bootstrap Normal 95% CI [0.00834,0.0199] [0.0083,0.0199] [0.00842,0.0199]
Bootstrap Percentile 95% CI [0.00847,0.0199] [0.00838,0.0199] [0.00847,0.0198]
BCa 95% CI [0.00845,0.0199] [0.00833,0.0199] [0.00854,0.0199]
Bootstrap Percentile-t 95% CI [0.00839, 0.0199] [0.00846, 0.0199] [0.00855, 0.0199]
  • comparison, delta s.e.: 0.00293, 95% CI: [0.00836,0.0199]

Experience level that maximizes log wages.

  • Application: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Average Effect of Experience

  • Experience level that maximizes log wages.

  • Bootstrap Hypothesis Test

  • Other Bootstrap Variants

Implement: Exp Level Max. Wages

  • We now implement bootstrap s.e. and CIs for experience level that maximizes log wage.

    • \(\theta = \left | \frac{\beta_{1}}{2 \cdot \beta_2} \right|,\)
    • \(\widehat \theta = \left | \frac{\widehat \beta_{1}}{2 \cdot \widehat \beta_2} \right|,\)
  • Key difference from before:

    • parameter involves ratio of coefficients,
    • estimator divides by estimated coefficient,
    • cause for worry?

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot2, R = B, stype = "i")


Bootstrap Statistics :
    original   bias    std. error
t1* 30.23942 4.915459    217.3218

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b, conf = 0.95, type = c("norm", 
    "perc", "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   (-449.04,  498.47 )   (  25.14,   52.12 )   (  25.22,   53.14 )  
Calculations and Intervals on Original Scale

Output from function boot

  • Note extreme skew to the right, extreme outliers.

  • explanation?

      R original bootBias bootSE bootMed
1 10000   30.239   5.5227 241.72  30.257
       V1           
 Min.   :    2.328  
 1st Qu.:   28.015  
 Median :   30.257  
 Mean   :   35.762  
 3rd Qu.:   33.556  
 Max.   :23085.192  

Limiting X-axis:

Results Using Delta Method

\[\sqrt{N}\left(\widehat \beta - \beta\right) \stackrel{d}{\rightarrow} N(0,\Sigma).\] Let \(\widehat \theta= |\frac{\widehat \beta_1}{2 \widehat \beta_2} |\), and \(\theta_0 = |\frac{\beta_1}{2 \beta_2} |\). Using delta method, we have \[\sqrt{N}\left(\widehat \theta - \theta_0\right) \stackrel{d}{\rightarrow} N(0,V)\] with \[V= a^{\prime}\Sigma a = \frac{\sigma^2_2}{4 \beta_2^2} + \frac{\sigma^2_3 \beta_1^2 }{4 \beta_2^4} - \frac{\sigma_{23} \beta_1}{2 \beta_2^3} \] where \(a=[0, -1/2\beta_2, 2\beta_1/\beta_2^2]\), \(\sigma_{j}^2\) is the \((j,j)\) element of \(\Sigma\) and \(\sigma_{jk}\) the \((j,k)\) element of \(\Sigma\).

Notice powers of \(\beta_2\) in the denominator of \(V\).

Results Using Delta Method

reg.1

Call:
lm(formula = I(log(wage)) ~ experience + I(experience^2), data = df)

Coefficients:
    (Intercept)       experience  I(experience^2)  
      1.5829802        0.0352680       -0.0005831  
a <- c(0,-1/(2*coef(reg.1)[[3]]),coef(reg.1)[[2]] /(2*(coef(reg.1)[[3]])^2))
delta.se <- as.numeric(sqrt(t(a)%*%vcovHC(reg.1)%*%a))
print(paste0("Delta method s.e.:  ",signif(delta.se,3)))
[1] "Delta method s.e.:  3.9"
print(paste0("Resulting Normal 95% CI: ","[",signif(maxexp-qnorm(.975)*delta.se,3),",",signif(maxexp+qnorm(.975)*delta.se,3),"]"))
[1] "Resulting Normal 95% CI: [22.6,37.9]"

Bootstrap Percentile-t CI

f.boot2.t <- function(d,i){ 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  est.b <- abs(coef(reg.b)[2]/(2*coef(reg.b)[3]))  
  a <- c(0,-1/(2*coef(reg.1)[[3]]),coef(reg.1)[[2]] /(2*(coef(reg.1)[[3]])^2))
  var.b <- as.numeric(t(a)%*%vcovHC(reg.1)%*%a)
  c(est.b,var.b)
  }
boot.out2.t<- boot(data=df, statistic=f.boot2.t, R=B, stype="i")

Bootstrap Percentile-t CI

boot.out2.t<- boot(data=df, statistic=f.boot2.t, R=B, stype="i")
boot.ci2.t <- boot.ci(boot.out2.t, conf=0.95, type=c("norm","perc","bca","stud"))
boot.ci2.t
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out2.t, conf = 0.95, type = c("norm", 
    "perc", "bca", "stud"))

Intervals : 
Level      Normal            Studentized     
95%   (-715.04,  760.31 )   (   9.19,   35.43 )  

Level     Percentile            BCa          
95%   (25.23, 50.42 )   (25.29, 51.55 )  
Calculations and Intervals on Original Scale

Repeated Simulations

boot.out.b2<- boot(data=df, statistic=f.boot2.t, R=B, stype="i")
plot(boot.out.b2)

summary(boot.out.b2)[1,]
      R original bootBias bootSE bootMed
1 10000   30.239   6.1939 218.67  30.243
summary(boot.out.b2$t)
       V1                  V2        
 Min.   :    6.459   Min.   : 6.667  
 1st Qu.:   28.023   1st Qu.:12.811  
 Median :   30.243   Median :14.736  
 Mean   :   36.433   Mean   :15.211  
 3rd Qu.:   33.470   3rd Qu.:17.120  
 Max.   :17844.393   Max.   :39.388  
boot.out.b3 <- boot(data=df, statistic=f.boot2.t, R=B, stype="i")
plot(boot.out.b3)

summary(boot.out.b3)[1,]
      R original bootBias bootSE bootMed
1 10000   30.239   5.5605 263.12  30.249
summary(boot.out.b3$t)
       V1                  V2       
 Min.   :    0.807   Min.   : 7.14  
 1st Qu.:   28.061   1st Qu.:12.84  
 Median :   30.249   Median :14.77  
 Mean   :   35.800   Mean   :15.16  
 3rd Qu.:   33.410   3rd Qu.:17.02  
 Max.   :25872.485   Max.   :34.31  

Repeated Simulations

boot.ci(boot.out.b2, conf=0.95, type=c("norm","perc","bca","stud")) # Simulation 2
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b2, conf = 0.95, type = c("norm", 
    "perc", "bca", "stud"))

Intervals : 
Level      Normal            Studentized     
95%   (-404.54,  452.63 )   (   9.32,   35.39 )  

Level     Percentile            BCa          
95%   (25.17, 50.54 )   (25.22, 51.13 )  
Calculations and Intervals on Original Scale
boot.ci(boot.out.b3, conf=0.95, type=c("norm","perc","bca","stud")) # Simulation 3
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b3, conf = 0.95, type = c("norm", 
    "perc", "bca", "stud"))

Intervals : 
Level      Normal            Studentized     
95%   (-491.02,  540.38 )   (   7.72,   35.48 )  

Level     Percentile            BCa          
95%   (25.13, 51.65 )   (25.21, 52.96 )  
Calculations and Intervals on Original Scale

Summary Repeated Simulations

Simulation 1 Simulation 2 Simulation 3
Estimated bias 5.5227 6.194 5.56
Bootstrap s.e. 241.72 155 186.35
Bootstrap Normal 95% CI [-449,498.5] [-404.5,452.6] [-491,540.4]
Bootstrap Percentile 95% CI [25.14,52.12] [25.17,50.54] [25.13,51.65]
BCa 95% CI [25.22,53.14] [25.22,51.13] [25.21,52.96]
Bootstrap Percentile-t 95% CI [9.19, 35.4] [9.32, 35.4] [7.72, 35.5]
  • comparison, delta s.e.: 3.9, 95% CI: [ 22.6, 37.9 ].

  • what’s going on? what would you report? how summarize?

Take-away lessons from application

  • We have investigated bootstrap inference for two parameters:

    • Average effect of experience: \[\beta_1 + 2 \cdot \beta_2 \cdot \mathbb{E}(\mbox{exp}),\]
    • Experience level that maximizes log wage: \[\left | \frac{\beta_1}{2 \beta_2}\right |.\]
  • How to summarize results for each parameter? What results would to rely on?

  • How are bootstrap results different for the two parameters? what do you think are the causes of the differences? what are the take-away lessons?

Bootstrap Hypothesis Test

  • Application: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Average Effect of Experience

  • Experience level that maximizes log wages.

  • Bootstrap Hypothesis Test

  • Other Bootstrap Variants

Test Using Asymptotic Normality

Consider

  • \(H_0: ~ \theta = \theta_0,\)
  • \(H_1: ~ \theta \ne \theta_0\).

Test based on asymptotic normality: \[ T_N = \frac{\widehat \theta_N - \theta_0}{\mbox{s.e.}(\widehat \theta_N)}, ~~~~~~~~~~~ T_N \stackrel{d}{\rightarrow} N(0,1) ~~ \mbox{under} ~ H_0\]

Test Using Asymptotic Normality

  • Compute asymptotic p-value based on asymptotic normality: \[ 2 (1 - \Phi(| T_N | ))\]

  • To test at level \(\alpha\), reject if p-value less than \(\alpha\).

Example

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

  • Consider null of linearity:
  • \(H_0: ~ \beta_2 = 0,\)
  • \(H_1: ~ \beta_2 \ne 0\).

Example

Using studentized test statistic and asymptotic normality: \[ T_N = \frac{\widehat \beta_2}{\mbox{s.e.}(\widehat \beta_2)}, ~~~~~~~~~~~ T_N \stackrel{d}{\rightarrow} N(0,1) ~~ \mbox{under} ~ H_0\]

test.stat <- coef(reg.1)[3]/sqrt(vcovHC(reg.1)[3,3])
p.value1 <- 2 * (1- pnorm(abs(test.stat)))
[1] "Test Statistic: -2.89"
[1] "P-Value: 0.00391"

Bootstrap Hypothesis Test

Again consider

  • \(H_0: ~ \theta = \theta_0,\)
  • \(H_1: ~ \theta \ne \theta_0\).

Again define studentized test statistic as: \[ T_N = \frac{\widehat \theta_N - \theta_0}{\mbox{s.e.}(\widehat \theta_N)}\] computed on original sample.

Bootstrap Hypothesis Test

Again define studentized test statistic as: \[ T_N = \frac{\widehat \theta_N - \theta_0}{\mbox{s.e.}(\widehat \theta_N)}\] computed on original sample.

  • Use same studentized test statistic as before, but use bootstrap distribution of studentized test statistic instead of normal distribution.

Bootstrap Hypothesis Test

  • Given \(B\) bootstrap replications, compute p-value
    \[\frac{1}{B} \sum_{i=1}^B 1\{ |T^*(b) | > |T_N| \},\] where \(T_N\) is test statistic on original sample and \[T^*(b) = \frac{\widehat \theta^*(b) - \widehat \theta}{s(\widehat \theta^*(b))}.\]

  • To test at level \(\alpha\), reject if p-value less than \(\alpha\).

Example

  • \(H_0: ~ \beta_2 = 0,\)
  • \(H_1: ~ \beta_2 \ne 0\).

Define:

\[ T_N = \frac{\widehat \beta_2 }{\mbox{s.e.}(\widehat \beta_2)}\] \[ T^*(b) = \frac{\widehat \beta_2^*(b) - \widehat \beta_2}{s(\widehat \beta_2^*(b))} \]

Example

f.boot.p <- function(d,i){ 
  reg.b <- lm(I(log(wage))~ experience+ I(experience^2),data=d[i,])
  test.stat <- (coef(reg.b)[3]-  coef(reg.1)[3])/sqrt(vcovHC(reg.b)[3,3])
  test.stat
  }
boot.test<- boot(data=df, statistic=f.boot.p, R=B, stype="i")
p.value.boot <- mean(abs(boot.test$t) >= abs(test.stat))

Example

Use studentized test statistic but bootstrap distribution of studentized test statistic.

f.boot.p <- function(d,i){ 
  reg.b <- lm(I(log(wage))~ experience+ I(experience^2),data=d[i,])
  test.stat <- (coef(reg.b)[3]-  coef(reg.1)[3])/sqrt(vcovHC(reg.b)[3,3])
  test.stat
  }
boot.test<- boot(data=df, statistic=f.boot.p, R=B, stype="i")
p.value.boot <- mean(abs(boot.test$t) >= abs(test.stat))
[1] "P-Value: 0.005"

Comparing

Figure 1: Normal Approximation

Figure 2: Bootstrap Approximiation

Other Bootstrap Variants

  • Application: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Average Effect of Experience

  • Experience level that maximizes log wages.

  • Bootstrap Hypothesis Test

  • Other Bootstrap Variants

Paired vs Residual vs Wild Bootstrap

  • We have done nonparametric bootstrap for OLS, also called “pairs bootstrap” in regression context.

  • Other bootstrap procedures for regressions:

    • Residual bootstrap,

    • Wild bootstrap.

  • Tradeoffs?

Non-i.i.d data

  • We have done bootstrap for i.i.d. data.

  • Can modify bootstrap procedure for

    • clustered sampling
      • straight forward modification, resample clusters
    • time series
      • less straight forward (block bootstrap).